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Abstract 

The onset of cohective behavior in a population of globahy coupled oscillators 
with randomly distributed frequencies is studied for phase dynamical mod- 
els with arbitrary coupling. The population is described by a Fokker-Planck 
equation for the distribution of phases which includes the diffusive effect of 
noise in the oscillator frequencies. The bifurcation from the phase-incoherent 
state is analyzed using amplitude equations for the unstable modes with par- 
ticular attention to the dependence of the nonlinearly saturated mode |aoo| 
on the linear growth rate 7. In general we find |aoo| ~ Vlil + l^D) where 
D is the diffusion coefficient and / is the mode number of the unstable mode. 
The unusual (7 -|- I'^D) factor arises from a singularity in the cubic term of 
the amplitude equation. 
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Typeset using REVT^ 



The onset of collective oscillations is a multi-faceted phenomenon of interest in physics, 
chemistry, biology and most recently neuroscience. One important class of models 

describes a collection of dissipative (limit cycle) oscillators that have some weak mutual 
interaction. For sufficiently weak couplings, the basic form of the uncoupled cycle persists, 
and the fundamental effect of the interaction is to alter the frequencies and evolving phases 
of the oscillators. When the coupling between two oscillators, f{6j — 6i), is uniform for all 
oscillator pairs, then the evolution of the phases is given by 

K ^ 

^^ = ^^ + ]vE/(%-^«) + e.W (1) 

for i = 1, . . . , N. If the oscillators are uncoupled {K = 0), then the phases simply evolve 
according to the unperturbed frequencies uji whose distribution is described by a density 
g{uj) characterizing the population. This normalized density is taken to have zero mean 
(cJ = 0); this can always be achieved by changing variables, 6i —* 6i — uJt, if necessary. 
For both physical and mathematical reasons, it is interesting to include in (|1]) the effect of 
extrinsic white noise ^j(t) perturbing the deterministic phase dynamics; this perturbation is 
defined by the ensemble averages < C,i(t) >= and < ^j(s) ^j(t) >= 2D6ij6{s — t). 

The form of the coupling function /(0) depends on the description of the underlying 
limit cycles and their mutual interaction, and will vary from one setting to another. 
Mathematically, since the coupling is necessarily 27r-periodic, we describe the general 
form by its Fourier expansion 

oo 

/(0)= E fne^'^'^- (2) 
n=— oo 

The simplest nontrivial possibility arises when the coupling is dominated by a single Fourier 
component, and the theoretical literature is largely focussed on the case /(0) = sin0 since 
this describes a strictly attractive interaction between oscillators with different phases. Early 
work by Kuramoto showed that, in the absence of noise {S,i{t) = 0), there was a critical 
coupling strength Kc above which a population g{uj) would begin to show frequency locking 
and partial ordering in the phases 6j. This transition was analyzed by an order parameter 
R defined as the time average of R(t): 
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1 ^ 

i?(t)e^^W = A ^ e^e^W. (3) 
i=i 

For large N ii K < K^, then R ^ 0, and for K > K^, R scaled like R ^ {K - K^)^/"^. 
In the limit ^ oo, this numerical result was also obtained analytically from a self- 
consistent calculation of R. For a population of identical oscillators {giyj) = S{uj)) in the 
presence of noise, Kuramoto analyzed the system of stochastic equations for the phases via 
the resulting Fokker-Planck equation which described the phase distribution. In this theory, 
the solution with R = becomes unstable for K > and a new state with R ~ {K — K^Y^'^ 
emerges. The value for Kc depends on the noise strength D in this case. Subsequently, work 
has generalized the Fokker-Planck approach to treat the phase dynamics for /(0) = sin0 
in populations with nontrivial frequency distributions gi^uj). These studies also show a 
bifurcation to phase ordered states with the same scaling R ~ {K — K^)^/'^. 

For couplings more general than = sin0, the properties of (|l|) are not understood 
and this is an interesting subject for several reasons. First, the couplings that are derived 
when a reduction to phase dynamics is actually carried out can easily have a more com- 
plicated structure. Secondly, recent results by Dai do indicate that as the form of /(0) 
is modified, the nature of the scaling exponent R ~ {K — K^)^ can change from the value 
1/2. 1^ Thus different forms of /(0) will correspond to different univerality classes. Daido 
specifically considers (|I|) without noise and applies his "order function" formalism devised 
as a generalization of Kuramoto's self-consistent calculation of R in the limit oo. |^ 

His treatment assumes the transition is triggered by the /±i components of the coupling 
and analyzes the self-consistent equation perturbatively to leading nonlinear order. The na- 
ture of the solution depends on a certain expression "Im(/_2C')" where (7 is a complicated 
function of the Fourier components of the order function. If this expression vanishes then 
the solution scales in the usual manner with f3 = 1/2, but if this expression is non-zero then 
the solution scales with /3 = 1. Thus a coupling with /1/2 7^ is predicted to produce tran- 
sitions with weaker phase ordering near onset than the transitions associated with the sin 
coupling. A third motivation for analyzing the transitions in (|I]) due to couplings of general 
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form is the unusual character of the bifurcations found in the Fokker-Planck description of 
these phase-ordering transitions: In the absence of noise, the unstable modes correspond 
to eigenvalues emerging from a neutral continuum at onset. IpHllI] This same feature has 



been noted in instabilities in other systems such as coUisionless plasmas, [|T2],|T3| ideal shear 



flows P^^TB|, solitary waves |]T^,[TB|, and bubble clouds. It is known that in some of 
these systems the nonlinear interactions between the unstable modes and the continuum 
can be singular in the sense that the amplitude equations for the modes become singular as 
the eigenvalue approaches the continuum. [|l2],|l3|,|l6l In these cases the singularities serve to 
alter the "expected" scaling behavior of the unstable modes. This is known not to occur for 
(HI) when /(0) = sin0, but Daido's results suggest this conclusion may depend crucially 



on the form of the coupling. 

In this Letter, I analyze the bifurcation from the phase incoherent state with i? = 0, in 
the limit N oo, within the framework of the Fokker-Planck equation P- p!T| , pO[] 

dp djpv) ^ 9V ... 
dt 39 39^' ^ ' 

The density p{9,uj,t) is defined so that N g{u!)p{9,uj,t) d9 du describes the number of os- 
cillators with natural frequencies in [u, uo + duj] and phases in [9, 9 + d9]. Thus p{9, u, t) d9 
denotes the fraction of oscillators with natural frequency u and phase in [9, 9 + d9] and must 
satisfy the normalization J^^ d9p{9,u,t) = 1 when g{iu) 7^ 0. In the limit 00, the 

deterministic part of the phase velocity (||) is expressed as an integral over the population 

/'27r POO 

vi9,uj,t) = uj + K d9' du'f{9'-9)p{9\uj\t)g{uj'), (5) 

Jo J -00 

and the coupling /(0) is described by its Fourier expansion (0). 

Equations and provide a continuum description of the oscillator population for 
which issues of stability and bifurcation can be analyzed in some detail. In this limit, 
Kuramoto's order parameter is given by 

/'27r /-oo 

i?(t)e^^W=/ d9' duj' p{9',uj',t)g{u')e'', (6) 

Jo J-00 



and the incoherent state {R = 0) is described by the uniform distribution po = l/^^r; this 
distribution is an equihbrium for (||) since v{9,Lj,t) = u + K/q at p = po- By setting 
p{9, u, t) = Po + ri{9, u, t), the dynamics can be rewritten for rj: 

^^=Cr,+M{r,) (7) 

in terms of the hnear operator 

Cr^ ^ -(co + Kfo)^ + ^ fj d9' dco'f'i9' - 9) r^{9', u', t) g{u:') (8) 



and nonhnear terms 
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M{7]) =K d9' dcu'ri{9',uj',t) g{u') 

Jo J~oo 



r^{9,u;,t)f'{9'-9)-^{9,u;,t) 



■ (9) 



In (1^) - @, = df /d(f), and note that the normahzation of p imphes J^^ d9r]{9, uo, t) = 0. 

The hnear stabihty of po, the onset of hnear instabihty, and the subsequent nonhnear 
bifurcation have been previously analyzed in detail for the specific case /(0) = sin0. [p|-[TT[| 
For an arbitrary coupling /(0), the generalization of this stability theory is summarized 
here as a prerequisite for the bifurcation analysis. The operator C acts independently on 
each Fourier subspace exp(m^); consequently the spectrum can be described by analyzing 
£\E' = A\E' for functions '^{9,uj) = iplu) exp{in9) with n > 0. In general this spectrum has 
both eigenvalues and a continuous component: For each mode number n = 1,2, . . ., there 
is a line of continuous spectrum at Re A = —n'^D; in addition C has eigenvalues when the 
function 

Ar,{z)^l + K /: r du;——f^ — - (10) 

j-oo uj + A Jo — z — inD 

has roots. More precisely, if for a particular mode number n = I, one finds A; {zq) = 0, then 

A = —Uzq is an eigenvalue of £ and \E'(6',co') = ipi^uj) exp{il9) is the eigenvector for A with 

The adjoint operator {C'^A,B) = {A,CB), defined via the inner product {A, B) = 
Jq^ d9 J^^du A{9,uj)* B{9,uj), has a corresponding eigenfunction ^ = iJj{lu) exp{il9)/2TT 
where 



and = X"^ . The normalization in (|12D assumes the root zq under consideration is sim- 
ple; this assumption can be relaxed if necessary but is characteristic of the codimension one 
bifurcations. This adjoint eigenfunction satisfies (\^, \E') = 1 and defines the projection, 



7] — > (\E', rf) \E', from rj onto the \I' component of rj. 
The example of a Lorentzian population, 

e r 1 

9[^) = - 

TX 



(13) 



Lcu2 + e2. 

provides an instructive illustration. For given values of (/, e), the solutions to A; {zq) = 
are easily located and one finds that C has an eigenvalue 

X = -l[lD + e + K Imifi)] - tlK[Re{fi) + /o] (14) 

whenever Klm{fi) < — e. Since e > and K > 0, the occurrence of these modes requires 
a coupling such that Im(/;) < 0. For /(0) = sin0, this is only possible for / = 1, but 
in general the mode number is not constrained. These modes are linearly stable when 
ID+e+K lm{fi) > and become linearly unstable for K > where Kc = —{lD+e)/lm{fi). 
For /(0) = sin0 and I = 1 this reduces to the familiar result = 2(e + D). |^ 

For D > 0, the resulting bifurcation for K > K^, can be analyzed by a center manifold 
reduction which yields an amplitude equation describing the time-asymptotic behavior of 
the unstable mode. We introduce this amplitude by writing t] in terms of the critical linear 
modes and the remaining degrees of freedom S: rj{0, u, t) = [a{t)'^{9, u) + c.c] + S{6, u, t) 
where (\E', S) = 0. In terms of (a. A) the evolution equation (|^) becomes 

a = Xa + {^,Ar{r])) (15) 
— = £5 + Afiv) - [(^, Xiv)) ^ + c.c] ; (16) 

The center manifold theorem asserts that any new solutions created by the bifurcation can 
be found by assuming S' is a function of a and a*: 
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oo 



S{9,u,t) = H{9,u,a{t),a{ty) = ^ Hn{uj,a{t),a{ty) e' 



(17) 



For these solutions r]'^ we have ri'^{9,uj,t) = [a{t)'^{9,u) + c.c] + H{9,u,a{t),a{t)*), and 
their dynamics is described by the two-dimensional flow 



The center manifold dynamics (|T^) and the function H in ([T7|) are both constrained 
by the symmetries of the problem. The group 0(2) is generated by rotations (3 ■ {9,uj) = 
{9 + [3,uj) and reflections k ■ {9,uj) = —{9,uj) which act on functions rj{9,uj) in the usual 
way: for any transformation x ^ 0(2), (x ' v){^j^) = vix^^ ' i^,^))- The operators £ and 
JV commute with rotations for arbitrary choices of g{uj) and /(0); in addition if g{uj) = 
g{—uj) and f{(j)) = —f {—(/)), then C and A/" commute with the reflection k. In the latter 
circumstance the bifurcation problem has 0(2) symmetry, otherwise the rotational symmetry 
alone corresponds to SO (2). 

The Fourier coefficients of H for n > are zero unless n is a multiple of /, the mode 
number of the instability. Rotational symmetry implies that the non-zero coefficients have 
the form Hi = aahi{uj,a) for n = / and if„(u;,a,a*) = a"^ hn{uj,a) 6n,mi for m = 2, . . . 
where a = |ap denotes the basic S0(2) invariant. The functions hn{uj, a) are unconstrained 
by the rotations, but must satisfy /i„(— cu, a)* = hn{uj,(r) when reflection symmetry holds. 
Similarly, rotational symmetry implies the amplitude equation (|1^) must have the general 
form a = p{cr)a where p{a) is a real-valued function if the reflection symmetry holds; 
otherwise p{a) is generically complex-valued. 

For K near Kc, we expand p{a) = po +pi(7 + . . . and hn{uj, cr) = hnfl{uj) + /i„_i(u;)cr + . . . 
and seek the leading (and presumably dominant) nonlinear terms in the amplitude equation 



a = Aa + (^,Ar(r/^)). 



(18) 



m 



a = a[\ + + ■■■]. 



(19) 



The calculation of pi from ([T8|) yields 
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Pi 



-27iiKl 



fi < 'ip, h2i,o > +f2i < -ip, ■ip* >< 9, h2i,o > 



(20) 



where the brackets denote an integration over u: < A, B >= du A* B, and the function 



h2ifi is determined self-consistently to be ||2T 



h2i,o{^) 



2l{uj + Kfo -zo- i2lD) 



du' 



2IK21 {zo) J -00 uj' + Kfo - zq- i2lD 



(21) 



For the special case /(0) = sin0, the instabihty arises for / = 1 and the f2i component is 
zero. Then the second terms in (|20| ) and ( |21| ) are absent and the results of Ref. [|1T| are 
recovered. 

When f2i 7^ 0, the new terms change the appearance of the bifurcation significantly due 
to the factor Kip^ip* > in (|20|). From (|Tl|) and (|T2D one sees that the integrand has poles 
uj± = {Q/l — KJq) ± i{l'^D + 7)// above and below the contour along the real axis; here 
zq = {Q + i'~f)/l is the root and 7 the linear growth rate from A = —Uzq. The 7 — > 0^ limit of 
< ip^ip* > in pi produces a pinching singularity when D is small or zero, and this singularity 
contributes a factor (7 + t^D) to the denominator when the integral is evaluated: 

//Jm(/0 



< ilj.ip >-- 



(22) 



The remaining integrals < ip, /i2«,o > and < g, h2ifl > in ([20|) are well behaved as 7 ^ 
since all poles lie in the same half-plane. 

The effect of the singularity (7 + PD) in ( p2D is clarified by scaling the amplitude 



a{t) = y^7(7 + PD) r(7t)e-*«(*) 

so that the equations for r(r) and C,{t) from (\[% are nonsingular: 

d/T r 

= r(r) 1 + Re (pi) (7 + PD) + ■ ■ 
dr L 



(23) 



d^ 
dt 



Q-lm (pi) 7(7 + + . 



(24) 
(25) 



Here t = 'yt is the slow time scale determined by the linear instability, and the coefficients 
in (|2^) and (pSD are now finite as 7 ^ O"*" even when D = 0. Assuming that as t ^ 00 the 
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instability saturates with the mode amphtude tending to a non-zero hmit r(r) Too, then 



the magnitude of this mode |aoo| = y 7(7 + /^-D) Too determines the scahng exhibited by 
the entrained state. For (7 + l^D) sufficiently small, the amplitude equation ( |2^ becomes 
independent of 7 and D and the scaling behavior of the entrained state follows the explicit 
dependence shown in (^3]). For D > 0, there is a crossover from |aoo| ~ 7 for 7 > l^D 
to |aoo| ~ y/l for 7 < l^D. In the noise- free limit, this crossover does not occur and the 
|«oo| ~ 7 scaling persists as the true asymptotic behavior. Since 7 ~ (K — Kc) near onset, 
these results determine the scaling exponent | 

Oool ~ (-^ — Kc)^ • III the special case D — 
and / = 1, these conclusions support Daido's findings: if /2 = then the singularity is 
suppressed and the (3 = 1/2 scaling occurs; if /2 7^ 0, then the / = 1 instability leads to an 
entrained state characterized by /5 = 1. 

Several important and related questions remain: When f2i 7^ and the cubic singularity 
occurs what are the singularities of the higher order terms in the amplitude equation? Does 
the amplitude scaling in ( pB]) suffice to control the higher order singularities if they occur? 
Finally, when f2i = are there higher order singularities that alter (3 from the value 1/2 
predicted by the cubic analysis? These issues will be addressed in a subsequent paper. p2 
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